/* boost_sobol.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
 */

#ifndef BOOST_RANDOM_SOBOL_HPP
#define BOOST_RANDOM_SOBOL_HPP

#include <boost/lambda/bind.hpp>
#include <list>
#include <bitset>
#include <vector>
#include <src/data/poly_mod_2.h>
#include <src/data/TUple.h>
#include <src/miscellaneous/OS_basic.h>
#include <src/algorithm/Eratosthenes.h>
#include <src/algorithm/random/Generator_Facilities.h>

namespace boost {
namespace random {


/*! \brief Sobol Generator.

	This is a Sobol generator. You are guaranteed that each Sobol generator on the heap is independant of the 
	Sobol generators still in memory.
*/	
template<OS_size MaxBit = 32,typename bitset = std::bitset<MaxBit>, class First_Algorithm = Erathostenes<poly_mod_2<bitset > > >
class sobol_generator : public Generator_facilities<sobol_generator<MaxBit,bitset,First_Algorithm>>
{
public :
  typedef OS_double result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
  static const OS_bool has_fixed_range = true;
#else
  BOOST_STATIC_CONSTANT(OS_bool, has_fixed_range = false);
#endif

	typedef typename First_Algorithm::Body							Body;
	typedef typename TUple<Body,MaxBit>								direction_number_type;
	class sobol_generator_base;
	typedef typename smart_ptr<sobol_generator_base>				smart_sobol_ptr;

private : 
class sobol_generator_base
{
public :
	typedef typename sobol_generator::direction_number_type			direction_number_type;
	typedef typename sobol_generator::Body							Body;

	sobol_generator_base(direction_number_type iv) : iv_(iv) { 
	seq_number = 0;
	ix=0;
  }
  sobol_generator_base() { 
	seq_number = 0;
	ix=0;
  }
  virtual ~sobol_generator_base() { 
		std::list<direction_number_type>& breakpoint = singleton_mgr::instance().unused_direction;
		breakpoint.push_front(iv_);
  }

  OS_double operator()()
  {
		unsigned long im = seq_number++;
		OS_int j=0;
		for (;j<MaxBit;j++) { 
			if ( !(im & 1)) break;
			im >>= 1;
		}
		ix ^= iv_[j];
		return ( ix.to_ulong()  ) * singleton_mgr::instance().fac;
  };

  direction_number_type iv_;
  OS_int seq_number;
  bitset ix;
};
/*! \brief Sobol Generator Manager.

	This class handles the creation of all Sobol generators.
	This is a singleton.
*/	
class sobol_generator_mgr
{
	typedef typename sobol_generator_base::direction_number_type		direction_number_type;
	typedef typename sobol_generator_base::Body							Body;
public :
	sobol_generator_mgr() : first_sobol((*this)()) {
		fac = 1.0;
		for (OS_int i = 0; i< MaxBit;++i)
			fac /= 2;
	}
	smart_sobol_ptr operator() ()
	{
		direction_number_type iv;
		if (unused_direction.size() == 0)
		{
			Body first;
			do
				first = ++first_;
			while (first < Body( (unsigned long) 3));
			OS_int degre = first.degree();
			iv[0] = 1;
			for (OS_int j=1;j<degre;j++) {
				unsigned long test = (1ul << (j+1) )* first_sobol();
				iv[j] = test-(test+1)%2;
			}
			for (OS_int j=degre;j<MaxBit;++j) {
				Body i = iv[j-degre];
				Body ipp = first;
				size_t deg =degre ;
				while(deg) {
					if (ipp[0] ) 
						i ^= (iv[j-deg] << deg);
					deg--;
				}
				iv[j]=i;
			}
			for (OS_int j=0;j<MaxBit;++j) {
				iv[j] <<= (MaxBit-j-1);}
		}
		else {
			iv = * unused_direction.begin();
			unused_direction.pop_front();
		}
		return smart_sobol_ptr( new sobol_generator_base(iv));
	};
	First_Algorithm					first_;
	std::list<direction_number_type> unused_direction;
	sobol_generator					first_sobol;
	OS_double fac;
	friend class sobol_generator_base;
};



public :
// Singleton over the Sobol Manager.
	typedef Singleton<sobol_generator_mgr > singleton_mgr;
	sobol_generator() {
		sobol_generator_ptr = singleton_mgr::instance()();
	}
	sobol_generator(const sobol_generator & rhs) : sobol_generator_ptr(rhs.sobol_generator_ptr) {};
	sobol_generator(smart_sobol_ptr & rhs) : sobol_generator_ptr(rhs) {};
	inline OS_double operator()() {
	  return (*sobol_generator_ptr)();
  }

  result_type min() const
  { return 0; }
  result_type max() const 
  { return 1; }
#ifdef BOOST_NO_OPERATORS_IN_NAMESPACE
    
  // Use a member function; Streamable concept not supported.
  OS_bool operator==(const sobol_generator& rhs) const
  { return sobol_generator_ptr == rhs.sobol_generator_ptr; }
  OS_bool operator!=(const sobol_generator& rhs) const
  { return !(*this == rhs); }

#else 
  friend OS_bool operator==(const sobol_generator& x,
                         const sobol_generator& y)
  { return x.sobol_generator_ptr == y.sobol_generator_ptr; }
  friend OS_bool operator!=(const sobol_generator& x,
                         const sobol_generator& y)
  { return !(x == y); }
#endif

protected : 
	smart_sobol_ptr sobol_generator_ptr;
};




} // namespace random

} // namespace boost

#endif // BOOST_RANDOM_SOBOL_HPP

